{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--BOOK_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"figures/PHydro-cover-small.png\">\n",
    "*This is the Jupyter notebook version of the [Python in Hydrology](http://www.greenteapress.com/pythonhydro/pythonhydro.html) by Sat Kumar Tomer.*\n",
    "*Source code is available at [code.google.com](https://code.google.com/archive/p/python-in-hydrology/source).*\n",
    "\n",
    "*The book is available under the [GNU Free Documentation License](http://www.gnu.org/copyleft/fdl.html). If you have comments, corrections or suggestions, please send email to satkumartomer@gmail.com.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Measure of Statistical Dependence](05.06-Measure-of-Statistical-Dependence.ipynb) | [Contents](Index.ipynb) | [Polynomial Regression](05.08-Polynomial-Regression.ipynb) >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.7 线性回归\n",
    "\n",
    "线性回归是一种用线性函数来模拟两个变量之间关系的方法。我们将使用`st.linregress`函数进行线性回归。我们首先使用已知的线性模型生成一些虚拟数据，并使用正态分布随机变量添加一些噪音。除了模型系数之外，`linregress`还提供相关性、p值和标准误差的估计。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.84221079098 7.09731123177\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "# 生成数据\n",
    "n = 100 #数据长度\n",
    "x = np.random.rand(n)\n",
    "y = 3 + 7*x + np.random.randn(n)\n",
    "# 执行线性回归\n",
    "b,a,r,p,e = st.linregress(x,y)\n",
    "print(a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们使用线性模型($y=3+7x+\\epsilon$)生成数据,而线性回归为($y=2.84+7.10x$)。拟合模型和真实模型的差异是因为噪音。当你增加更多的噪音时，你会看到拟合模型更加偏离现实。图5.13显示了真实线($y=3+7x$),损坏测量($y=3+7x+\\epsilon$),拟合线($y=2.84+7.1x$),和拟合线的预测区间。拟合线和真实线匹配合理。预测区间也相当合理。\n",
    "\n",
    "预测的$Y_pred$的方差由下式给出:\n",
    "\n",
    "<center>$\\sigma^2 = E[(Y_pred-\\widehat{Y})^2] = \\sigma_{\\epsilon}^2(1+\\frac{1}{n}+\\frac{(X_0-\\overline{X})^2}{\\sum_{i=1}^{n}(X - \\overline{X})^2})\\qquad$(5.1)<center>\n",
    "    \n",
    "其中，$\\sigma_{\\epsilon}^2$是由剩余方差的经典无偏估计的$s^2$估计的。$\\sigma_{pred}^2$是使用t分布与n-2自由度产生的预测区间(因为$s^2$是一个估计量)。在$Y_pred$置信区间由下式给出,\n",
    "\n",
    "<center>$PI = \\sigma_{pred}*z\\qquad$(5.2)</center>\n",
    "\n",
    "其中，$PI$是预测间距，$z$是t分布在$\\alpha$显著水平上的值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = y - a - b*x # 拟合数据和测量数据的误差\n",
    "x1 = np.linspace(0, 1) # x轴绘制PI\n",
    "# 拟合误差的方差\n",
    "e_pi = np.var(eps)*(1+1.0/n + (x1-x.mean())**2/np.sum((x-x.mean())**2))\n",
    "# z值使用t分布和dof = n-2\n",
    "z = st.t.ppf(0.95, n-2)\n",
    "# 预测间距\n",
    "pi = np.sqrt(e_pi)*z\n",
    "zl = st.t.ppf(0.10, n-2) # z在0.1\n",
    "zu = st.t.ppf(0.90, n-2) # z在0.9\n",
    "ll = a + b*x1 + np.sqrt(e_pi)*zl # 10 %\n",
    "ul = a + b*x1 + np.sqrt(e_pi)*zu # 90 %"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最终，我们可以绘制出真实线、拟合线、带噪音的损坏测量和预测区间。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4VGX2wPHvm05CEAhFEJKANGnSpUiTRV3E3kVARFhR7OsPFV11XVbXzrquCqIgxQoii9gVQZoEEUSQ3nsNhBSSzPn9cRMIYSaZmczMnXI+zzNPZu7cmfveBO65bzuvERGUUkpFrii7C6CUUspeGgiUUirCaSBQSqkIp4FAKaUinAYCpZSKcBoIlFIqwmkgUEqpCKeBQCmlIpwGAqWUinAxdhfAHTVq1JD09HS7i6GUUiFl2bJlB0SkZnn7hUQgSE9PJyMjw+5iKKVUSDHGbHVnP781DRlj3jHG7DPGrCqx7QVjzB/GmJXGmE+NMVX9dXyllFLu8WcfwUTg0lLbvgFaikhrYB3wqB+Pr5RSyg1+CwQiMg84VGrb1yJSUPRyMVDPX8dXSinlHjv7CG4HPvT2w/n5+ezYsYPc3FwfFimyJSQkUK9ePWJjY+0uilIqgGwJBMaY0UABMLWMfYYDwwFSU1PPeH/Hjh0kJyeTnp6OMcZfRY0YIsLBgwfZsWMHDRo0sLs4SqkACvg8AmPMYKA/MEDKWBVHRMaJSAcR6VCz5pmjn3Jzc0lJSdEg4CPGGFJSUrSGpULX1KmQng5RUdbPqS7vM1UpAa0RGGMuBUYBPUUk2wffV/FCqZP096lC1tSpMHw4ZBddVrZutV4DDBhgX7lChD+Hj74PLAKaGmN2GGOGAv8BkoFvjDG/GmPe9NfxlVIRZPToU0GgWHa2tV2Vy5+jhm4WkToiEisi9URkgog0EpH6ItKm6HGnv46vytarVy+dpKfCx7Ztnm33lTBpjoqcXENh8gcrS0FBQfk7KRWOnAwoKXO7LxQ3R23dCiKnmqNC8NoSGYHAT3+wLVu20KxZM+644w5atmzJgAED+Pbbb+nWrRuNGzfm559/5vjx49x+++107NiRtm3b8tlnn538bPfu3WnXrh3t2rVj4cKFAOzevZsePXrQpk0bWrZsyfz58wGoXLnyyeN+8skn3HbbbQDcdtttPPjgg/Tu3ZtRo0a5PF5OTg433XQTrVu35sYbbyQnJ6dC565UUBkzBhITT9+WmGht95dwao4SkaB/tG/fXkpbvXr1GdtcSksTsULA6Y+0NPe/w4nNmzdLdHS0rFy5UgoLC6Vdu3YyZMgQcTgcMnPmTLnyyivl0UcflcmTJ4uIyOHDh6Vx48aSlZUlx48fl5ycHBERWbdunRSf44svvij/+Mc/RESkoKBAjh49KiIiSUlJJ4/78ccfy+DBg0VEZPDgwXLZZZdJQUGBiIjL47300ksyZMgQERFZsWKFREdHy9KlS884J49+r0oFkylTrP/Txlg/p0zx7/GMcX5dMca/x/UAkCFuXGNDIulchfmx/bBBgwa0atUKgBYtWtCnTx+MMbRq1YotW7awY8cOZs2axYsvvghYw163bdtG3bp1GTlyJL/++ivR0dGsW7cOgI4dO3L77beTn5/PVVddRZs2bcotw/XXX090dDQAX3/9tdPjzZs3j3vvvReA1q1b07p16wqfu1JBZcCAwI4QSk21WhecbQ8xkREI/PgHi4+PP/k8Kirq5OuoqCgKCgqIjo5m+vTpNG3a9LTPPfXUU9SuXZsVK1bgcDhISEgAoEePHsybN4/PP/+cgQMH8vDDDzNo0KDThnaWHuuflJR08rmIOD0e6PBQpXxqzJjTh6yC/5uj/CQy+gjsaD8scskll/Daa68hRXPnli9fDkBmZiZ16tQhKiqKyZMnU1hYCMDWrVupVasWw4YNY+jQofzyyy8A1K5dmzVr1uBwOPj00089Pl6PHj2YWtQnsmrVKlauXOmfE1YqUgwYAOPGQVoaGGP9HDcuJOctREYgsPEP9sQTT5Cfn0/r1q1p2bIlTzzxBAB33XUXkyZNonPnzqxbt+7kXf3cuXNp06YNbdu2Zfr06dx3330APPfcc/Tv35+LLrqIOnXqeHy8ESNGkJWVRevWrXn++efp1KmTn89cqQgwYABs2QIOh/UzBIMAgCm+cwxmHTp0kNJj3tesWcN5551nU4nCl/5elQofxphlItKhvP0io0aglFLKJQ0ESikVrApOBOQwGgiUUiqYOBynnu9eEZBDaiBQSqlgkHsUFr4G/z4fDm22ttVsEpBDR8Y8AqWUClZHd8HiN2DZRMg7CundIb9obkLCWQEpggYCpZSyS94x+E9H68Lf/Croeg+c0y7gxdCmoQq4/fbbqVWrFi1btjy57dChQ/Tt25fGjRvTt29fDh8+DMD06dNp0aIF3bt35+DBgwBs3LiRm266yZayK6VsIAKb58O3T1uv45Oh/6tw73K4/l1bggBoIKiQ2267jS+//PK0bc899xx9+vRh/fr19OnTh+eeew6Al156icWLFzNo0CCmTZsGwOOPP84zzzwT8HIrpQKssABWzYDxvWFSf1g+GY4fsN5rfT1US7e1eBoIKqBHjx5Ur179tG2fffYZgwcPBmDw4MHMnDkTsHIP5eXlkZ2dTWxsLPPnz6dOnTo0btw44OVWSgXQ7pXwWjv4ZIjVIdz/Vbj/N0iqYXfJTgqbPoIb31p0xrb+reswsEs6OScKue3dn894/7r29bi+Q30OHT/BiCnLTnvvw7908aoce/fuPZkCok6dOuzbtw+AJ598kksuuYS6desyZcoUbrjhBj744AOvjqGUCnJZ+6xO4LptoHoDSGkEl/wTmvazFscKMmETCIJd37596du3LwCTJk2iX79+rF27lhdffJFq1aoxduxYEksnxlNKhZYDG2DRa/Dr+9bFf8QCqx9g4Ay7S1amsAkEZd3BV4qLLvP96klxXtcASqtduza7d++mTp067N69m1q1ap32fnZ2NpMmTeKrr77i4osv5rPPPmPatGlMnTqVYcOG+aQMSqkA2/Ur/Pg8rJ0D0XHQ5mboMtJKchkCgq+OEuKuuOIKJk2aBFh3/ldeeeVp7z///PPcd999xMbGkpOTgzGGqKgosksveaeUCm4OBxTkWc8PrINtC6HHX+GBVXD5WKgROv1/YVMjsMPNN9/M3LlzOXDgAPXq1ePpp5/mkUce4YYbbmDChAmkpqby8ccfn9x/165dZGRk8NRTTwHw0EMP0blzZ6pWrXqyU1kpFeTyc2DF+7DodWh7K1z4ALS4GppdBnFJ5X8+CGkaanUa/b2qsDB1qrWI/LZt1kqEY8a4XivA3X2zD8HSt2HJW5B9AOq2hZ6PQNNL/XsuFeBuGmqtESilwsvUqacvIbl1q/UazrzAe7LvzLtg3RfQ+GLoei+kXxgyfQDl0T4CpVR4GT369HWEwXo9erRn++5cBh8PgSPbre0XPQ53LYYBH0OD7mETBEBrBEqpcLNtm/vbnW1rHANd98P4iyD+LDj/JqhaH85ueea+YUIDgVIqvKSmWk08zraXtW8UMCwJzo6GLAMX/wPaD7bmAYQ5bRpSSoWXMWOg9ORMY6BfvzP3/ftoaF800scB/FEAswuhzb+h60jfBIGpUyE93ZpRnJ5uvQ4yGgiUUuFlwAAYPPj0NnwRmDTp1EX4yDb48lHY/iT0j4aW9az9t5wND7wFtw7yTVmKO6O3brXKUNwZHWTBQANBBY0dO5aWLVvSokULXn31VUBTUStluzlzrAtvSdnZMOZR+GQojG0DP4+D8/rDnT/Bb9utCWJbtrgeZuoNTzqubeS3QGCMeccYs88Ys6rEturGmG+MMeuLflbz1/EDYdWqVYwfP56ff/6ZFStWMHv2bNavX6+pqJWyW+lO4MSi2sHmHbB1AXQeAfetgGvGwdmtAleO8rbbxJ81golA6ZkWjwDfiUhj4Lui1yFrzZo1dO7cmcTERGJiYujZsyeffvqppqJWym6pqdbVrXUs3JkEA4v6DGqnWimgLxkDZ9ULTDk82W4Tv40aEpF5xpj0UpuvBHoVPZ8EzAVG+eSA71525rYWV0GnYXAiG6Zef+b7bW6BtgPg+EH4qFSb4JDPyz1ky5YtGT16NAcPHqRSpUrMmTOHDh06aCpqpeyUmwl/7QWbZ0CygX2FsOgEJCVaHcnRsYEry5gxp09YA6sje8yYwJXBDYEePlpbRHYDiMhuY0yt8j4QzM477zxGjRpF3759qVy5Mueffz4xMa5/pZqKWqkAWDMbDnwKtZrCrF2waLd1B/5WGWkm/KX4eO6mu7CLiPjtAaQDq0q8PlLq/cNlfHY4kAFkpKamSmmrV68+Y5vdHn30UXn99delSZMmsmvXLhER2bVrlzRp0uS0/Y4fPy69e/eWEydOSK9evSQzM1PeeOMNGTdunB3FPk0w/l6VKtPu30SmDxdZ/Jb1Oj9XZOcv9pYpSAAZ4sa1OtCjhvYaY+oAFP3c52pHERknIh1EpEPNmjUDVkBPFTf7bNu2jRkzZnDzzTdrKmql/E0ENv4Ak6+BN7vBmv/BiWPWezHxVkI45bZANw3NAgYDzxX9/CzAx/e5a6+9loMHDxIbG8vrr79OtWrVNBW1Uv42+wFY9i5Urg19/gYdbodKIT0I0VZ+S0NtjHkfq2O4BrAXeBKYCXwEpALbgOtF5FB536VpqANHf6/KrzxJD11S3jH45T1oeS0knw1bfoJDm6D1jVYNQDllexpqEbnZxVt9/HVMpVQQ8yTlc7Gju2HJm5DxLuRlWhf9jndYKaDTLwxMuSOAzixWSgWGq1m2t956Zg4ehwM+GwmvtoKF/4Zze8Ed31tBwJUQyOkTrEI6EPirWStS6e9T+VVZs2mLawfv/Mt6HRUF4oAOQ+CeZXDDe1CvvevPe5LTRwPGGUJ2qcrNmzeTnJxMSkoKJowWiLCLiHDw4EGOHTtGgwYN7C6OCkfp6c7TQxugeQx0jYe60XDXEqjVzDffnZZm5Q8qVrp5CqwJXuPGBd/Yfh9wt48gZANBfn4+O3bsIDc316ZShZ+EhATq1atHbGwAZ16qyFH6IhwDtI+DznFQNQoOFMLifFiUCbEJnn13VNSZSebAyijqcJx67W7ACBO2dxb7W2xsrN65KhUM3B0JVLzt8cdgyzaINdAnHnYVwhe5sK7AuiB7GgTA/cVoQiQJXKCFdB+BUspmnrTN7/sDKi+EJ86DyZPBVIL/ZMHEbCsIVCQHj7PFaJx9X4gkgQs0DQRKKe+Vl29fBDbPh6k3wH8vgN+mQ+0WcON1Vrt8tVSr+SYtrWLt9AMGWJ9PSyv7+9wNGBEmZPsIlFJBoLy2+VUz4JMhkFgDOg23hn8mpQS+nCV5O6ktBLnbR6A1AqWU90o3qcQCneKg59nW66Z/hsv/DQ+sgl6j7A8CYF30t2zxz4pkzoTAcNWQ7SxWSgWB4nz7UTlWAOgQB5UMJLew3o+tBO0H21tGO3kzm9oGWiNQSnlvwAD41/VwfzJcGAd746Dew/DQN3aXzD88vbsPkTWLtUaglPKMiLXub81mkFQD+t4E66tDl7sh5Vy7S+c/Zd3dg/N+hxAZrqqBQCnlnsICWPMZLHwNdi2H3o9Dz4eh6aXWI9y5uru/7z7IyXEeINyd32AzbRpSSpVvyTh4rS18cjvkHoX+r0DXkXaXKrBc3cUfPOi6+SdEhqtqIFBKOZebeer55h8huQ7cOBVGZlgLwcRW8v0xg3mEjad38du2uT+/wWY6j0Apdbp9f8Ci1+C3T+DOBVCjEZzIhrjE8j9bEcGeEM5V+SpVsmoFpQVB/iKdR6CUcp+zGcBtbz111+/vIAD2jbBxtxbi6u5+7NiQaP4pi3YWK6Ug+xBMuRbik6HXY/bMALZjhI2n4/wHDHBdOwnh2craNKRUJMrLguVTYMdSuG6CtW3LT3BOe/+0/bvDjhTRYZ6WWpuGlFJnOrYHvn0aXmkOX46CzB3WKCCw1gC2KwiAPSNsgnyc/96jgVlvRQOBUpFi4w/wSkv46RVo0BOGfgNDv4KEKt59X0VG+Dj7rB0jbIIsLbXDIcxbtx+Hw2qpOZCVF5gDi0jQP9q3by9KKQ85HCKbfhRZ9431Oi9L5ItHRQ5sqPh3T5kikpgoYnUzW4/ERGu7Pz/ra0FSluy8Apm8aIv0fuEHSRs1W77/Y69PvhfIEDeusVojUCrcFBZYQz/H9YJJl8P8F63tcUlw6T99kwbC1QifwYPLryEEU/4dV7UQCMh8huwTBbz89Vq6Pvcdj89cReWEGP59c1subFTDL8dzRTuLlQonv38KX/8NMrdBSiPoMhLOv8n3bf/GlL+PqzkA7q4vbBdn8wWMgTvvhG7dfDI6KDMnn7MqxVJQ6KDXi3NpdnYVhnVvQKcG1THu/G7dFPaL1yulihzdBbGJUKmqtRDM0reh6z3Q+BLrousPMTFQWFj+fs5G3wT7SB1X5QOIjYX8/FOvPZjwJiL8tOEA4+dv5o/dR5k/qjfxMdFknyggMc4/I/l11JBS4W7v7/DpCHi1Nfw83trW4moYMsdaEMZfQQDcCwLgfPRNsOffKWvEUMkgAG41aeUVFPJxxnb+PHY+Ayf8zOpdRxnUJY3Cog5hfwUBT2ggUCrUbJoLk6+BN7rC6plW3p9W11nv+bBZwani0T7uSk09c4QQBHf+HW9yCpVh2dbDPPzJSkTg+etas+CR3oy8qHFQBIBiwVMSpZRrDsepO/xF/4U9v8FFj0OHoZBY3ffHc7auL5zZdl6WxETo18/5zN1x44KjGciZMWNg4EDn/RjOlAocWw8eZ8JPm0lOiOHhS5rRpWEKHw7v7PP2f1/SPgKlglnuUfhlEvw8Dgb/D6qlw9HdUKkaxCb455ieJlcD666+Xz+YM+f04DF6dHD3B7hy113w5punB4O4OOu1kz4CueUWlm09zPj5m/h69V5iogwDLkjjqStaBL7sJWhnsVKhLHMnLHkDlk2CvKOQ3h0ufQ7Obun/Y5fVWepMWaN9gn2EUFlc1YqcjBp69dt1vPrtes6qFMutnVMZ3CWdWlX8FKg9ENSBwBjzAHAHIMBvwBARcTmXWgOBiih5x+DFplCQA82vskYAndMucMd3dfF2pay7+2AZIeTsol6BPomsvAI+XLqdLg1TaF63Cuv2HmPxpoNc175eULX9uxsIAl5iY8w5wL1AcxHJMcZ8BNwETAx0WZQKCiKw8Xtr8Ze+f7cygF4+Fup3gmppgS+Pq+UVU1JOX5IRyh/tM2aM82amQI4Q8jTDaBl2Z+YwccEWpv28jWO5Bdz/p8Y0r1uFJrWTaVI72ccFDxy7Rg3FAJWMMTFAIrDLpnIoZZ+CE/Dr+/BGN5hyDaz4EI4XtcG3vt6eIACuh3eOHev5aJ9gWKHLRzOZn5i5iu7/+oG3f9pMzyY1mXl3N+7/UxMfFtRG7uSh8PUDuA/IAvYDU13sMxzIADJSU1N9kndDqaAwZYpIu3NEHqgs8mQVkTHNRJZPFcnPs6csaWkixlg/i3PsuNoeiow5PZdQ8cOYMj9WWOiQ+ev2i8PhEBGRl79eK3//3++y7eDxQJTaJ3Az15AdTUPVgCuBBsAR4GNjzK0iMqXkfiIyDhgHVh9BoMuplM8d2QbT34N7n4WCbGhcCWblwO5tkCbQJi6w5SmvySRYxvVXlKumLhfzBXLzC5nxy04m/LSJjfuPM3loJ7o3rskDfcPk7t8JO5qG/gRsFpH9IpIPzAC62lAOpQJj16/wyVAY2wZ+fdm68J4APsyBjYX2JVwLpuRv/uTmTObsEwW88s06uj33PY99+hsJsdG8emMbOjcM8EptNrAjEGwDOhtjEo01u6IPsMaGcijlX1sXwsT+MK4nrPsKOo+AyZnO9y1vIZSK5P53JcgXZfGZcvopsvIKAIiJiuKjjO2cX78q04ZdwOx7LuSqtucQGx3+CRgC3jQkIkuMMZ8AvwAFwHKKmoCUCnn5uSAOa7H3I9vh4EZrJFD72yDhLKg6FTLdb6YAfDrq5YxjetBkEtJKNXWJCIs2HuDtogRwP/5fb+Jiovj6gR4kJ8TaWFB72BLqRORJEWkmIi1FZKCIBGgZHhVR/HEX7Ur2IZj3ArzaypoFDNDyWrh/JXS7zwoC4F3CNX814QR78jc/yC90MHP5Tvq/9hO3jF/Ciu1HuKFjffILrcltkRgEQHMNqXDlr7vo0g5tgsVvWAvB52dDoz9Z4/8Bop389yo+tieTm/zVhONNWULc4k0Huf/DXzm3ZhLPXtOKq9ueQ0JstN3Fsp2mmFDhKVAzWidfA5vnQesboMvdUNsPuWWCZXaup3w8m9cb2w9l8+6CLZxVKZb7/tQYEWHBhoN0PTeFqKjgTADnS7oegYps/riLdhTCmtlWB3DmDmvbpc/B/b/BVf/1TxCAijXhBLJ5rPRxhw+3ApjIqRpZgI6/YvsR7p72Cz1f+IH3Fm05uQi8MYYLG9eIiCDgCQ0EKjy56vD0piP0RDYsnQD/6QgfDoDDW60HQM0mUKWO9+V0h7ezc+28GNs4NHXst+u58vUFzFu7n2HdGzJ/VG+eOb4CatSwfn/GWM8DFRRDgTuzzux+tG/fvqIT7FSkmTJFJDHx9JmkiYmez5A9kS3yQhNrBvBbvUR+my5SkO+fMvtaWprzGbVpaf4/tpezed1WYuZzdsNG8t7Yj2TdnqMiIrJq5xEZP2+jHM05cWrf2NgzyxIXF9ozpt2AmzOLbb/Iu/PQQKC84kmahJL7tqon8trtp95b+LrI5p9EilINhAx/X4zL4s8gVBTk9yVWlRe73ypt7pkqaaNmy8svfeJZWQIVFG3kbiDQzmKliptQauRBlzhoGmPNcGnyHAy+y+7Sec9XnczedPq6WtzGFwnn0tN5utHFTG3Tj/zoaP60fgnDls6kY3QWxtl5lZVWOxTWRagAdzuLbb/bd+ehNQLlV63riQxNspp/Hq4s0iteJMmE/t2iL5rHKvIdPkxc53A4ZOnmg1YCOGPkuR6DZXTfEbKxWt3yazpaIyj3Gmv7Rd6dhwYC5XO5x0T2rbWeJxiRYUkiHWJFYgLchOJvFb0Y29nPICJ5+YUyfdl2+fOr8yRt1GxZuOGASEqK8zKlpDj/Eu0jKPcaqxPKVGQ5tgeWvAUZE6BqKvxlPtROhfFhmmqhollEbcpHlHOikIkLtzBx4Wb2Hs2jca3KPH9ta9qmVvX8y4rP/777Tq25nJJira8QxpPnPKHDR1Vk2L8WZt4Fr7SEBa9Cw15w2StWG7FdqRbsGuPvCV8Ow3VDbn4hYP1ZJvy0mUa1KjNxSEe+fqAHN3Ssb80CPnTI+YddbQfrgn/gwKn6wIEDGgRKcqfaYPdDm4aUVxyOU0M9l70n8o+zRT7/q8jBjWfuG+iFWHw1vNXfXJVzxAif/r5+2XpI7pqyTHo8/73kFxSKiMihLBcL9djcXBVK0FFDKmIV5sOqGbDwNWhzC3S5Cwry4MRxSKxud+ksoZQ2ovSooX79YNKkCo8IKnQI367Zy9vzN7F0y2GSE2K45YJU7r2oMUnxZbRa+3NEUphxd9SQBgIVPnIzYdlEWPwmHNsFNZtB78eg+ZV2l+xMroY0hsJwRh8Fsblr93Hbu0s5p2olbr+wATd2rE/lsgJASUGQxygUaCBQkWfq9bD+a2jQA7reC+f2sS64wciOGoGvLp5eBrH9x/J4b9EWKsfH8Jee5+IoqhFc1KwWMRGw+Isd3A0EOmpIha5dy2HRf62FX6rUse7+e4+Gum3sLlnZpk6FrKwzt/uzg9qXabk9XNBm/d5jvD1/M58u30m+w8F17eoBEBVluLjF2Z4dW/mFhmEVWhwOa9nHif1hXC9Y+wXs+c16r27b4AoCzkYFFV+Qi4cxFktJqXgbd1mjkHyZBK5fP7e3v/7DBvq+Mo/PVuzkho71+O7Bnrxw/fmeH1P5ldYIVHBy1ozhyIeMB6CaA44bSLsG7nj11OpfwcTVHXilSmdekAEqV654ECjrjt+X8wHmzHG5Pb/QweyVuzi/XlUa1qzMhY1qUOgQbu2cRvWkOM+PpQKi3D4CY8xIYKqIHA5Mkc6kfQQRpuRFrZKBRtGwtqhd+sIoOFAIvxdAQhCPFHHVB+BKRTuJy+tzqGifRMnA7OSakRmfxPttLmXiFSPYczSXkb0b8ddLmnp4EsrXfNZZbIz5B3AT1mLz7wBfSYB7mDUQRJj0dDi6DbrEQ5tYq976ahYcdfLPLhiHW0LZic6cqeh5lNeBW5Ehl84+W8IL3Qcysf3lHI9PpEvDFIb3aEjPJjV18Zcg4LMVykTkcaAxMAG4DVhvjPmnMebcCpdSqdIOb4XO++CeytA2Flblw3+POw8C4PdUB15zNfM2JcU/s5jLmwHs7eI24LR/YW2NNIr/IsfjKtF38zJmNzrG+8M707tZLQ0CIcatzuKiGsCeokcBUA34xBjzvB/LpiKFoxCO7rKexydDvTiYf8KqBczKhQNlNJkEaz4gV2krxo71/oLszfFKBpgBA6xah8Nh/XT3mEXB1oHh23M7ccPNz3LJ0NfJOKc5GMOTG77i1etb0fKOmyp2Dso25XYWG2PuBQYDB4C3gYdFJN8YEwWsB/7Pv0VUYSsvC36dCoteh8QUGPa9NfO3zWvwxp2QXaIWEBdnNX3k55/aFoh8QN4qvsi6Grfv636N8o5XAXnpDfmkSmMmdLySTSn1OSdzH49//zbNEgrB4UDv/UOfO30EfwcmiMgZPU3GmPNEZI2/CldM+wjCTMkMoLmZUP8C6DISzrvcuksG56OGQGeTBlChQ4iOMuRMnkrXjCjqHdnLHUs/pd/aBcQmxAdvR706SWcWq+AjYl3oM96Bzx+CZv2h6z1Qv5PdJXNfBKQ22Lg/i7fnb2b5tsN8fm93oqMMu96dRp2nH8OE8XmHI51ZrIKDCGyaayWAa/pn6DQMzr/ZSgNdvaHNhfOQL2fn+qo8PgpKIsKSzYd4e/4mvl2zj7iYKK5tV4/sEwUkJ8RSd8gtMOQWH5+AChYaCJR/FJyAVdNh0X9g7ypIqgXNr7Dei60UekEAyp6dG8hAMHXq6YusQIVHuQinAAAb8UlEQVSD0rz1Bxj8zs9UT4rj3j6NGdQljRqV431UYBXstGlI+ccHA+CP2VDzPOg6ElpdDzEhfmEJhoyh5Yzpd3c+QlZeAR8u3U5ctGFgl3QKHcKny3fSv3Uda/EXFRZ8No9AKbcc3gJfPgpZ+6zXXe6GAdPhrkXQ9tbQDwLgv9W6PFmpzFmtpKRy5lXszszh2Tlr6PLsdzwzezULN1q1iugow3Xt62kQiFDaNKQqZkeG1f6/ZhaYKEjtYjUBpXW1u2S+N2aM89m5FRnC6mm/Q3kT6MoISu8u2MyYz9fgEOHPreowrHtD2tT3Yg1gFXZsqREYY6oaYz4xxvxhjFljjOliRzlUBRScgHf7wdt9YOMPVv7/+3871Q8QjioyO9cVT7OCllX7KBWURIS5a/ex/ZD1/a3OOYuB1XL5cdbfeP3WDrTp3iY410lWgefOepa+fgCTgDuKnscBVcvaX9csDhJ5x0XWfX3q9ed/FVn0X5Hco/aVKdQZ43z9XWOc7+9sDWEQSUk5uW5wbn6BfLh0m/R9ea6kjZot/5yz2vVnvVknOdDrOyuvEaxrFhtjqgArgIbi5sG1s9hmWfvg5/Gw9G3IOQz3rYBqaXaXKjy4mxW05FDR6kXrLh86dMaw0fHzNjFu/ib2H8uj2dnJDOvekMvPr0tcTJRvVkXT9YJDSjB3FjcE9gPvGmOWG2PeNsYkld7JGDPcGJNhjMnYv39/4EsZajzpcHTX0d0w6x54pSXMewFSO8Ntn0PVIM3vEww8/Tu4kyOo+OK7dat1H3/wIOTkwOTJsGULe/pfe3LXjfuzOK9OFaYMvYAv7uvOte3rWUEAfLMmgS8XuFHBw51qgy8fQAesxHUXFL0eCzxT1me0aagcvqryi4g4HCLZh6znx/aJPJsqMus+kf3rfVvmcOTt36G8ppa0NKfNR8va9ZQ7J2dIg0dmyy9brb9ZfkGh6+O4+B5JS3P/HD1tylK2Ioibhs4GFotIetHr7sAjInKZq89o01A5fFHlL8yH3z+1JoDFVIKhX1nb83OsCWCqfP5akL7E/IVCE8U3jS7g7U5XkVGvBVUSYhjQOY0h3dKplZxQ9vf4olnHX+eo/CJom4ZEZA+w3RhTvHxRH2B1oMsRVipS5c/NhAX/hrHnw4xh1oW/zc2nJkiFWhDwRxOZu3y5HGRJqaknc//nxsTxf/3uY0/lFJ5c9jGLHu3DqEublR8EwDejntxpylKhx51qg68fQBsgA1gJzASqlbW/Ng2VoyJV/iXjRJ6sIvLuZSJ/fCFSWEbTQrDzZRNZ8fd5MjrGF00vpew7misvvfSJXH/r81KI1SzzR400yU+qbN9oHR01FDJws2nIlkDg6UMDQTk8uQDuWCby8RCRZZOs13lZIjuXB7a8/uLLC7E3QcWHgWj93qMy6pMV0nj0HEl/ZLbcMeZTOdKomV58lUfcDQSaayhclJWJ0uGA9V9ZM4C3LoC4ZOj9qJUGIpz4MheQt23hPsgIumjjQW4ev5j4mCiua1+PoRc2oGHNyh59h1Kg6xGokj4eAr/PgCr1oPOd0G4wJFSxu1Rl8+aC6suOzAAmmCsodDBn1R4KCh1c064eBYUOJvy0meva1yNFM4CqCtD1CCJZ1n5r9a9Ow62lH9sNhGaXQfMrITrW7tKVz9u8/77MBZSa6jyo+HCN5OIMoO/8tJmdR3K4oEF1rmlXj5joKP7S81yfHUep8mj20XCyfx3MuhdeaQFzn4WN31vbz70IWl3nOgjYOdLGGW8nLfkyF5CfR8d8nLH9ZAbQc6pWYvygDrw/rLNPvlspT2mNIBwUnICPBsG6LyAmAdrcYrX/12hc/meDbdUtqNgwzAEDfFNuPywGv2b3UaonxVG7SgLnVK1EjyY1NQOoCgraRxCqCvNhx9JT6Z4/HWHl/+l4ByTVcP97gnGCUDCWyUsiwvz1Bxg/fxPz1x9gWPcGjL6sud3FUhFC+wjCVW4mLJsES96EY7utBHBVU+HqN7z7Pn9NgqoIf+T9t8HM5Tt588eN/LHnGLWS4/m/S5syoJMm61PBR/sIQkXWfvjyMXi5BXzzhLXm703vWyOBKqKiq26527/gST+EP/L+B8jxvIKTz+et248IvHBda+aP6s1dvRpxVmIIdNaryOPOZAO7HxE9oSwvy/qZuUvkH2eLfDJUZOcvvvv+ikyCKu+zxTNQi5OS+WrGbxDadvC4PD3rd2n+xBfy244jIiJyLDdfHA6HzSVTkQw3J5RpjSAYORzwx+fwzp9h2o3Wtip14ME1kNsbul7tuxE+Fbn7Lmt0T8nUyXDmmPwwSV28cscRRk77hZ4v/MB7i7ZwSYuzqRxvtbhWjo/BGGNvAZVygwYCf/J0WOaJbFg6Af7TAT64BTK3Q5NLT01gmvH56Xnpi0f4+CIYbNli5bcHGDjQvfKW1b9Q3iLrZX0+RGSfKGDA+CX8uHY/w7o3ZP6o3rx8YxvSa5yxvIZSwc2daoPdj5BsGvKmyWXRG1YCuLd6iqz8WKQg//T3K5JLp7xEYd6Ut6zyuMpb76NkbHbIOVEg05ZslTsnZ5xs8lm08YAczTlhc8mUcg5NOmczdy7ae9eIzLxbZMWH1gW3UapIWoxIWqrzC7C3i4K4c5F3p7ylg8mIEa6/19X3uRNkgiy75cGsPBn77Tpp/8zXkjZqtvQbO0/2H8u1tUxKuUMDgd3Kumhv/EFk8rXW3f8ztUT+Pci9u3FvawTufK68IOMqmIwY4fyi7Wz/4mOUdXH3dSrpCvpl6yFp+vgcSRs1W257Z4ksWL9fO4BVyNBAYDdXF99ba1gB4PlzReb+SyTrgPsXeG8vku7UJMorg6v3o6PLbm7y5M5+yhTr+2xuRsrYcki++X2PiIjk5RfKU7NWydo9RwN2fKV8RQOB3Yov2gmIdI0TiS+6aL81yloL4ETOqX09afLxptnE3WafsoKMO23+FV0EpvTxPWn+qqCCQod88dtuuea/CyRt1Gy55JUf9c5fhTwNBHY7tFnklStEHqti1QB61XF9kfTDylancbcmUVaQKa/Nv6JlLu/7/Vgj+Hb1Hun5/PeSNmq2XHj3RHm3/eWSdW4T2/smlKooDQR2yc8T+XCQyFNVRZ5OEZnxF5FdK8r+TCDaxSvaAVveHXtF79zLqnH4oY9g/7FcOZiVJyIiP/yxV654aqbMbn2R5Jsovx5XqUDSQBBIhQUi2zNOvf5wkMjXfxPJ3On+dwTZSBmnSpbRVVt+dLR3ZS+rD8KHv4sN+47JI9NXSuPRc+Sfc1aLiIjD4RCHv2tlStnA3UCg2UcrIi8Lfp0Gi1+HI9vh/pVwVgVz/4SK0umrS0pM9Dw3kLPv8+Z7XFi65RBv/biJb9fsJa7EEpDnFi8BGcAVyZQKFM0+6k/HD8Ki1yDjXcg9AvU6wp+ehuQ6dpcscIovzoMHQ2Hh6e8Vp4/w5ALuh/z/DocQFWWleHhv0VaWbT3EvX0aM6hLGjVKLwEZgBXJlApWWiPwRH4uxCZYd/+vtbPSP3S9B+p38s33+2Dh84ALwjvp7BMFfLR0O+8s2ML4QR1oenYy+47lkhwfS6W4aOcf8nONRCk7aI3AV0Rgw7ew6D9gomHgDKhaHx78A5JSfHecYFwpzB1BdCe971gu7y3cyuTFW8nMyad9WjVy863aSq3khLI/7IcaiVKhQpPOuZKfC7+8B//tDFOvg/1roUH3U3e/xUHAV+v9ertOr938vLavu/IKCrn4lXm8PncDXc9NYfqIrkwf0ZXzPVkGsjj5nsNh/dQgoCKE1ghcyZgAXz0GtVvB1W9Bi2sgJu70fXx5Fx+MK4W5w6Y7aRFh8aZDfL16D3/r35z4mGjGXNWKFnWraPZPpTykfQTFDmywRv806AEtroacI7D7V2jQ02rvdsaXa+uG0Tq9/lRQ6GDOqj2Mn7eJ33ZmkpIUx//uuZC6VSvZXTSlgo67fQSR3TQkAlt+gmk3WWsALJ8KhzZb71WqCg17wbRprpt+fHkXP2YMxJWqccTFhdw6vf60bu8xer4wl3vfX87xvAL+eXUrFjxykQYBpSoospuGZgyH3z6CxBTo+X/Q8Q6oXOvU++U1/fi6o7R07SwEamv+tvdoLtsPZdMhvTqp1RNpXrcKT13Rgj7Nap0cGqqUqiB3Zp3Z/fDZzOKcIyILXhPJLcokufp/IksniJzIdr5/ebNNK7reb8mZxCkpFZ/ZGgqzk930x+6j8tBHv0qjxz6XHs9/7zoBXBids1K+hpszi22rERhjooEMYKeI9PfrwQ5vhSVvWqOATmRZ6/+2vBbOK+ew5TX9eNtR6qym4WkZ3PnOUBh+WsqK7Ud4+Zt1/LhuP5Vio7mlUypDL2zofO3fMDlnpexmW2exMeZBoANQpbxA4HVncUEefPoXWP0ZYKDlNdDlbqjb1r3P+6sD19X3OhMdbQ1nLC/IhHBnc36hg/xCB4lxMXyzei+PzviN27qmMeCCNKolxbn+YAifs1KB4G5nsS2BwBhTD5gEjAEe9FsgAHj/Fkg5Fy74i+d5gPw129TVbNzylHXsIJzhW56svAI++Hkb7/y0mWvb1+Ohi5vicAgnCh0kxLqYAVxSCJ6zUoEU7DOLXwX+D0j2+5Funub9Z/01Rt5VJ3PlypCTc2bunmJl5fAJohm+5dmTmcu7Czczbck2juUWcEGD6nRMrw5AVJQhIcqNIAAhdc5KBbOADx81xvQH9onIsnL2G26MyTDGZOzfvz9ApXPCH7NNnc3GjYuDvDzXQaCYqz6DIJnh645nZq9m/LxN9GxSk8/u7saHf+lCjyY1Pf+iEDpnpYKaOz3KvnwAzwI7gC3AHiAbmFLWZ4J+PQJvuDtqyJNRREE4gsbhcMj8dftl0IQlsmHfMRER2bw/S7YdPO6bAwThOSsVLAiF9QiMMb2Av4o/+whChTv9BiGUDTO/0MHnK3czbt4mVu8+Ss3keJ6/rjW9m9Yq/8NKKZ8I9j4CVZqr9m53Rw0FkYJCB5e8Mo9NB47TqFZlnr+2NVe2rUt8jJtt/0qpgLI1EIjIXGCunWUIGmPGhHQ+/D2ZuXyxajdDujUgJjqKQV3SSE1JpFcTnQGsVLCL7FxD4Ls00hU1YIB10U9Ls4Y/pqVVPAgE4NzW7D7Kgx/9yoX/+p5nZq9m0/4sAG7r1oCLmtXWIKBUCIjs7KPBvCpVRVcr8/O57TySwyPTVzJ//QES46K5sWN9bu/WgPrVE8v/sFIqIIJ6Qpmn/BYIgnVmqi8u4n44t/xCBzsP55BeI4mcE4Vc9foCrmhTlwEXpFI1sYwZwEopW2ggcEdZI3XS0uxbstAXF3Efzro9lpvPBz9v550Fm0mIjebbB3sSHWWsYWeu1mpQStlO1yNwh6sZqMZYF2KRU4nMKtq+7kl7vS/WOXB1bh7Mut2Tmcuzc9bQ9dnvGTNnDWkpifytf3OKm/01CCgVHiI7EDibmWrMmXfSFV07uLipx93g4oOLeEVm3Toc1vn/su0w4+dvomfTmswa2Y0Phnehd7NaGgCUCjfuzDqz++HXmcWlZ6a6mtFrjPfHKG9dA2dl8nadg7LOrYzPF88AHjhhibz6zToRESkodPhuBrBSKuAIhZnF7grozGJ/dCB7015f0VFDbio9A7hG5Xju7dOIQV3SfX4spVRgaR+Bt7xtUimrD8Cbph5/JLtz4rEZv3H/h79yotDBv65txYJHemsQUCrCaIqJ0rxJPV3eSlmuZg3bkCWzOAX0LZ1SSUtJYnDXdC5teTa9m+oMYKUilTYN+YI7zUkBaupx5Y89Rxk/bzOzVuyk0CE8d01rbuhYP2DHV0oFns4j8JY3F+wgXilLRBg+eRnfrN5LpdhobuhQj6EXNiQ1RWcAKxXuNPuoN7xdDD3IVsoqKHSwYONBejapiTGGc2tW5vyLzyp/DWClVETSGkFJ3o4YCpKcRVl5BXy4dDvv/LSZnUdy+N/IC2lV76yAHV8pFVx01JA3XM3c3bq17NnA/sgcWpZSI5Qy35vGv778g67Pfsczs1dzTtVKjB/UgRZ1q/jn+EqpsKI1gpJc1QiKBUNm0hK1j+zYeBLz88g6qzrd755Il+Z1Gda9IW1Tq9lXPqVU0NDOYm84a+JxJi3NttXCJD2dxY4qjO90Ndur1uarCSOJQjh2blOSN/wR8PIopYKXdhZ7o+QcgrJqBu52IvtQQaGDL3/fw7ie97OyTmNSjh9h8C+zyY+OJr6wgORN6wJSDqVU+NE+gtKKZ/SmpZW9X0UT0bmjRF/A1z2vZeS05WQlVeGfX77Ggjdv596FHxBfWGDta9MIJaVU6NNA4IqzVBOleZIW2kP7J07jxXe+473qLUGEvov+x4T/Pce3ZHDL+vkkFJw4tbNNs5SVUuFBA4ErJUcCueKHu/AN+7J4ZPpKuv2eyOsdrmZtTev4sY5C+qz+iag5cwI7QkkpFfa0j6AsxRfXIUMgP//09+LifH4X/u/v1vPyN+uIj4ni+pXfMHTpTBoe3nX6Ttu2WeXSC79Syke0RlCe0aPPDAIAyckVvhgXOoQvftvNriM5AFzQoDr39mnMgkcuYszaz88MAlCxWognq6QppSKGBoLyuOoHOHTI66/MOVHI5EVbuOiluYyY+gsfZWwH4IKGKTzYtwk1KsdXaIUxpzxdJU0pFTE0EJTHF8tGlvCf79fT7V/f88Rnv1MtMY43BrTjnosan7mjr2crjx595vyIQIx8UkoFvcgMBJ40kfjgznxPZu7J55sOHKddalU++ksXPr2rK39uVYdoV+sA+HJxGlc1Gz+OfFJKhYbICwSeNpGUd2deRlBZtvUwd05eRtfnvmPVzkwAXrjufN4e3JFODaoHdhF4H9dslFLhI/JSTPhyTWInKSkciUl8++IExlGfjK2HOatSLLd2TmVItwZW279dgiRDqlIqcDTFhCu+bCJx0u6eVSA8uCGaqrVyefLy5tzQoT5J8UHwa/ZmCU6lVETQGkExb2oEUVEciU9iStt+LKnfkvc++hsGWF2rAU12bSAmOvJa3pRSwSNo1yMwxtQ3xvxgjFljjPndGHNfQAvgo2GZ2w9l89RVD9FlxERe7DGIaIeDo/FJADSv5NAgoJQKGXa0WRQAD4nIL8aYZGCZMeYbEVkdkKP7oIlk0caDDHh7MdFNe3LF6h8ZtvBjmh0oqmVo3h+lVIgJeCAQkd3A7qLnx4wxa4BzgMAEAvA4RYOIMHftfnLyC+nXqg7t06pxz0WNublTKmfPzoQVH8NBo+3uSqmQZGsfgTEmHZgHtBSRo672C9jCNKWcKHDw2a87GT9/E+v2ZtEutSoz7uoW8HIopZQ3gn7UkDGmMjAduN9ZEDDGDAeGA6TaMNb985W7+fvs39l7NI9mZyfz0vXnc/n5dQNeDqWU8jdbAoExJhYrCEwVkRnO9hGRccA4sGoEgSjX7swc4qKjSKkcT2J8NI1qVeb5686nR+MagZ38pZRSARTwQGCsK+oEYI2IvBzo4zuzZvdRxs/bxKwVuxjavQGP/vk8ejetRe+mtewumlJK+Z0dNYJuwEDgN2PMr0XbHhOROYEuyMINB3hz3ibmrdtPYlw0A7ukcesF5SxRqZRSYcaOUUM/Aba1szgcQlRRkrdpP29j9a6jPHxJU269II2zEmPtKpZSStkmCHIfBMbxvAI+WLqdd37azLtDOtKkdjJPXdGCyvExJMRG2108pZSyTdgHgn3Hcpm4YAtTFm/laG4BHdOrkZfvALA3CZxSSgWJsA4EufmFXPzKPDJz8rmk+dkM79mQdqnV7C6WUkoFlbAOBAmx0Yy5qhXN61ahQY0ku4ujlFJBKawDAcBlrevYXQSllApqmiJTKaUinAYCpZSKcBoIlFIqwmkgUEqpCKeBQCmlIpwGAqWUinAaCJRSKsJpIFBKqQhn61KV7jLG7Ae2evnxGsABHxYnFOg5RwY958hQkXNOE5Ga5e0UEoGgIowxGe6s2RlO9Jwjg55zZAjEOWvTkFJKRTgNBEopFeEiIRCMs7sANtBzjgx6zpHB7+cc9n0ESimlyhYJNQKllFJlCJtAYIy51Biz1hizwRjziJP3440xHxa9v8QYkx74UvqWG+f8oDFmtTFmpTHmO2NMmh3l9KXyzrnEftcZY8QYE9IjTNw5X2PMDUV/59+NMdMCXUZfc+Pfdaox5gdjzPKif9v97CinLxlj3jHG7DPGrHLxvjHG/Lvod7LSGNPOpwUQkZB/ANHARqAhEAesAJqX2ucu4M2i5zcBH9pd7gCcc28gsej5iEg456L9koF5wGKgg93l9vPfuDGwHKhW9LqW3eUOwDmPA0YUPW8ObLG73D447x5AO2CVi/f7AV8ABugMLPHl8cOlRtAJ2CAim0TkBPABcGWpfa4EJhU9/wToY4wxASyjr5V7ziLyg4hkF71cDNQLcBl9zZ2/M8AzwPNAbiAL5wfunO8w4HUROQwgIvsCXEZfc+ecBahS9PwsYFcAy+cXIjIPOFTGLlcC74llMVDVGOOz5RfDJRCcA2wv8XpH0Tan+4hIAZAJpASkdP7hzjmXNBTrjiKUlXvOxpi2QH0RmR3IgvmJO3/jJkATY8wCY8xiY8ylASudf7hzzk8BtxpjdgBzgHsCUzRbefr/3SPhsmaxszv70sOh3NknlLh9PsaYW4EOQE+/lsj/yjxnY0wU8ApwW6AK5Gfu/I1jsJqHemHV+OYbY1qKyBE/l81f3Dnnm4GJIvKSMaYLMLnonB3+L55t/Hr9CpcawQ6gfonX9TizunhyH2NMDFaVsqyqWLBz55wxxvwJGA1cISJ5ASqbv5R3zslAS2CuMWYLVlvqrBDuMHb33/VnIpIvIpuBtViBIVS5c85DgY8ARGQRkICVjyecufX/3VvhEgiWAo2NMQ2MMXFYncGzSu0zCxhc9Pw64Hsp6oUJUeWec1EzyVtYQSDU246hnHMWkUwRqSEi6SKSjtUvcoWIZNhT3Apz59/1TKxBARhjamA1FW0KaCl9y51z3gb0ATDGnIcVCPYHtJSBNwsYVDR6qDOQKSK7ffXlYdE0JCIFxpiRwFdYow7eEZHfjTF/BzJEZBYwAasKuQGrJnCTfSWuODfP+QWgMvBxUb/4NhG5wrZCV5Cb5xw23Dzfr4CLjTGrgULgYRE5aF+pK8bNc34IGG+MeQCreeS2EL+pwxjzPlbzXo2ivo8ngVgAEXkTqy+kH7AByAaG+PT4If77U0opVUHh0jSklFLKSxoIlFIqwmkgUEqpCKeBQCmlIpwGAqWUinAaCJRSKsJpIFBKqQingUApLxhjOhblhU8wxiQVrQXQ0u5yKeUNnVCmlJeMMf/ASm9QCdghIs/aXCSlvKKBQCkvFeXCWYq17kFXESm0uUhKeUWbhpTyXnWsXE7JWDUDpUKS1giU8pIxZhbWCloNgDoiMtLmIinllbDIPqpUoBljBgEFIjLNGBMNLDTGXCQi39tdNqU8pTUCpZSKcNpHoJRSEU4DgVJKRTgNBEopFeE0ECilVITTQKCUUhFOA4FSSkU4DQRKKRXhNBAopVSE+3+PNfaep2ZoAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x191daa2fb38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(x,y,'ro', label='measured')\n",
    "plt.plot(x1,ll,'--', label='10%')\n",
    "plt.plot(x1,ul,'--', label='90%')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
